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It is shown that quadrupole susceptibility can be detected in Gd compounds contrary to our textbook knowledge that 
Gd 3+ ion induces pure spin moment due to the Hund's rules in an LS coupling scheme. The ground-state multiplet 
of Gd' i+ is always characterized by 3=112, where J denotes total angular momentum, but in a j-j coupling scheme, 
one / electron in j=H2 octet carries quadrupole moment, while other six electrons fully occupy j=5/2 sextet, where j 
denotes one-electron total angular momentum. For realistic values of Coulomb interaction and spin-orbit coupling, the 
ground-state wavefunction is found to contain significant amount of the j-j coupling component. From the evaluation 
of quadrupole susceptibility in a simple mean-field approximation, we point out a possibility to detect the softening of 
elastic constant in Gd-based filled skutterudites. 

KEYWORDS: Quadrupole, Gadolinium, Filled skutterudites, Wigner-Eckart theorem. 
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1. Introduction 

Recently, multipole phenomena have attracted much atten- 
tion in the research field of strongly correlated /-electron sys- 
tems. 13 In general, multipole indicates spin-orbital complex 
degree of freedom, which is considered to be active in /- 
electron materials due to the strong spin-orbit coupling of / 
electrons. However, when orbital degeneracy is lifted, for in- 
stance, due to the effect of crystal structure with low sym- 
metry, only spin degree of freedom often remains. Thus, /- 
electron compounds with high symmetry are important for the 
research of multipole phenomena. In this sense, filled skut- 
terudite compounds LJ1T4X12 (Ln: lanthanide, T: transition 
metal, X: pnictogen) with cubic structure of Th point group 
have provided us ideal stages for multipole physics. 4 

Here we pick up Gd-based compounds. In general, Gd ion 
takes the trivalent state in compounds, leading to the 4/ 7 sys- 
tem. In an LS coupling scheme, due to the Hund's rules, 
electron configuration is determined so as to maximize to- 
tal spin angular momentum S, indicating 5=7/2 for Gd 3+ . 
Since seven / orbitals are singly occupied, total orbital angu- 
lar momentum L is zero in this case. Then, the total angular 
momentum J of the ground state is given by J=S=1I2. The 
ground-state multiplet is generally affected by a crystalline 
electric field (CEF) potential, but the octet of Gd ion is not 
changed at all, since the CEF potentials act only on charge, 
not on spin, and the octet with J=S=1I2 and L=0 is spher- 
ically symmetric state. Namely, in the LS coupling scheme, 
there is no CEF splitting in the ground-state multiplet and we 
cannot expect to observe the response of higher-rank multi- 
pole in Gd-based compounds. 

Note here that the LS coupling scheme is exact in the limit 
of U 3> A, where U denotes the magnitude of Hund's rule in- 
teraction and A is spin-orbit coupling. In actual materials, U 
is larger than A, but X/U is not equal to zero. Thus, the effect 
of a j-j coupling scheme should appear, more or less, in the 
ground-state wavefunction. For Gd 3+ ion, the ground state in 
the j-j coupling scheme is also characterized by J=7/2, but it 
is composed of fully occupied j=5/2 sextet and one electron 
in j=7/2 octet, where j indicates one-electron total angular 
momentum. Since one / electron in the j=7/2 octet possesses 
orbital degree of freedom, we can observe the CEF level split- 



ting even for Gd compounds in the j-j coupling scheme. For 
actual rare-earth compounds, U is about a few eV, while A is 
several thousand Kelvin, indicating that X/U is in the order of 
0. 1 . Namely, the actual materials are in the intermediate sit- 
uation between the LS and j-j coupling schemes. Thus, the 
ground state of Gd compounds includes the contribution of 
the j-j coupling scheme and the quadrupole susceptibility is 
expected to be significant in contrast to our knowledge on the 
basis of the LS coupling scheme. 

In this paper, it is shown that quadrupole susceptibility can 
be detected in Gd compounds, although it has been simply 
believed that Gd 3+ ion induces pure spin moment due to the 
Hund's rules. It is found that the j-j coupling component ap- 
pears in the ground-state wavefunction for realistic values of 
Coulomb interaction and spin-orbit coupling. With the use of 
the CEF parameters determined from PrOs4Sbi2, 5 ~ 7 we eval- 
uate the quadrupole susceptibility for Gd-based filled skut- 
terudite compounds and discuss the detection of the softening 
of elastic constant. 

The organization of this paper is as follows. In Sec. 2, we 
define the local /-electron Hamiltonian and explain the defi- 
nition of multipole in the form of one-body operator. In Sec. 3, 
we discuss the CEF states and determine the parameters from 
the experimental results of PrOs 4 Sb 12 . In Sec. 4, we evalu- 
ate quadrupole susceptibility in a simple mean-field theory 
and discuss the detection of the elastic constant for Gd-based 
filled skutterudites. Finally, in Sec. 5, we summarize this pa- 
per. Throughout this paper, we use such units as /i=/cb=1. 

2. Local Hamiltonian and Multipole Operator 

2.1 Local Hamiltonian 

First we define the local /-electron Hamiltonian as 

-ffioc = Hq + H so + Hqef, (1) 

where Hq denotes Coulomb interaction term, H so is a spin- 
orbit coupling term, and Hqef indicates crystalline electric 
field (CEF) potential term. Among them, Hq is given by 
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where a=+l (—1) for up (down) spin, fi ma is the annihilation 
operator for / electron with spin a and z-component m of 
angular momentum 1=3 at an atomic site i, and the Coulomb 
integral I mi m 2 , m3 m 4 is given by 

6 

^2F k Ck{m 1 ,m 4 )c k (m2,m 3 ). (3) 

fe=o 

Here F k indicates the Slater-Condon parameter and is 
the Gaunt coefficient. 8 Note that the sum is limited by the 
Wigner-Eckart theorem to fc=0, 2, 4, and 6. Although the 
Slater-Condon parameters should be determined for the ma- 
terial from the experimental results, here we simply assume 
the ratio among the Slater-Condon parameters as physically 
reasonable value, given by 

F° = IOC/, F 2 = 5U, F 4 = 3U, F 6 = U, (4) 



where U denotes the scale of Hund's rule interaction among 
/ orbitals, which should be in the order of a few eV. 
The spin-orbit coupling term H so is expressed as 
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(5) 



where A is the spin-orbit interaction, and the matrix element 
£ is given by 



Cm,(j:rn,(j — 771(7/2, 

(m+a,-a;m,a = \[W + 1) = m ( m + 

and zero for other cases. 

Finally, the CEF term Hqef is expressed as 



Hcef - ^2 B m . m ,fl ma f, 



(6) 



(7) 



where B m ^ m i denotes the CEF potential for / electrons from 
the ligand ions, which is determined from the table of Hutch- 
ings for angular momentum 1=3. 9 For the filled skutterudite 
compounds with Xh symmetry, 10 B mm i is expressed by us- 
ing three CEF parameters, B®, B§, £>|, as 

B 3 .3 = 5-3,-3 = 1S0B 4 + 180B 6, 

B2.2 = B_ 2 ,-2 = -420B*} - 10805°, 
B hl = = 60B2 + 2700Bg, 



B, 
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B2-2 
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300B^ 

5-3, -1 = 
5-2,0 = - 
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notation, 11 we define £??, _Bg, and B? as 



,. Following the traditional 
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Wx/F(A), 

W(l-\x\)/F(6), 

Wy/F\6), 



(9) 



where x and y specify the CEF scheme for Th point group, 
while W determines the energy scale for the CEF poten- 
tial. As for F(4), F(6), and F*(6), we choose F(4)=15, 
F(6)=180, and F* (6)=24 for ^=3. 9 



2.2 Definition ofmultipole operator 

Next we define multipole as spin-orbital density in the form 
of one-body operator from the viewpoint of multipole expan- 
sion of electron density in electromagnetism. On the basis 
of this definition of the multipole operator, we have devel- 
oped microscopic theories for multipole-related phenomena. 
For instance, octupole ordering in NpC>2 has been clarified by 
evaluating multipole interaction with the use of the standard 
perturbation method in terms of electron hopping. 1214 We 
have also discussed possible multipole states of filled skut- 
terudites by analyzing multipole susceptibility of a multior- 
bital Anderson model based on the j-j coupling scheme. 15-26 
Quite recently, a microscopic theory for multipole ordering 
from an itinerant picture has been developed on the basis of a 
seven orbital Hubbard model with spin-orbit coupling. 27 

The multipole operator T is expressed in the second- 
quantized form as 



T 



(k) 



q : Tn<7 : m' a' 

where k is a rank of multipole, q denotes an integer between 
— k and k, 7 is a label to express Oh irreducible representa- 

(k) 

tion, G-y.q is the transformation matrix between spherical and 

n(k,q) 
ma, 

Wigner-Eckart theorem as 2 



cubic harmonics, and T^£' m(jl can be calculated by using the 

„28 



x (j>|&ns!>0V|&n'«y). 



(11) 



Here 1=3, s=l/2, j=£±s, /1 denotes the z-component of j, 
(j fj,\f p! j" p,") indicates the Clebsch-Gordan coefficient, and 
1 1 1 j) is the reduced matrix element for spherical tensor 
operator, given by 



{j\\T (k) \\3) = 



l_ / (2j + fc + l)! 
2 fe \/ {2j-k)\ ' 



(12) 



We note that k<2j and the highest rank is 2j. Namely, for 
/ electrons, we can treat the multipoles up to rank 7 in the 
present definition. 

In this paper, we define the multipole operator in the one- 
body form on the analogy of the multipole expansion the- 
ory of electromagnetic potential, but some readers may cast 
questions about the validity of the present definition. Since 
the multipole operators have been traditionally defined by the 
Stevens' operator equivalent in the space of the Hund's-rule 
multiplet characterized by the total angular momentum J, 29 
the present definition seem to be different from such a tradi- 
tional one at a first glance. Concerning this point, we provide 
detailed comments in the following. Hereafter the Stevens' 
operator equivalent is simply called the Stevens operator. 

First of all, we reply to the above question in a simple man- 
ner. In one word, the present definition of multipole operator 
is related to the Stevens operator through the Wigner-Eckart 
theorem as 



(n,J,M\fW\n,J,M') 

= c^j(U,X)J2G^Ol(J;M,M% 



(13) 
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Fig. 1. (Color online) Coefficients c of the multipole operators for n=7 
obtained from the diagonalization of Hq + H bo . 



where n is the local /-electron number, J is the total an- 
gular momentum of the ground state multiplet, M is the z- 
component of J, |n, J, M) denotes the degenerate ground 
states of Hq + H so , c denotes a proportional coefficient de- 
pending on n, J, k, U, and A, and Oi (J; M, M') is the matrix 
element of Stevens operator Ol, given by 



O q k (J; M, M') 



(J\\TW\\J) 



(JM\JM'kq). 



(14) 



\/2J + 1 

Here the reduced matrix element is given by eq. ( 1 2), in which 
j is replaced with J. We can numerically calculate all the ma- 
trix elements of the Stevens operator for arbitral numbers of 
J and k, but for some low-order multipoles, it is convenient to 
use the symmetrized polynomials of J Q , 30 ' 31 where J a (a=x, 
y, and z) is the a component of total angular momentum J. 

In Fig. 1, we depict the coefficients c for the case of n=7, 
corresponding to Gd 3+ ion. We note that the coefficients do 
not depend only on the ratio of X/U, but even if the values 
of U and A have been individually changed so as to keep 
the ratio, we have not found significant difference in the re- 
sults. Note also that in the present case of n=7, the ground 
state multiplet is characterized by J=7/2 and the highest mul- 
tipole is rank 7 in any case. As naively expected, we find that 
only dipole moment remains in the limit of the LS coupling 
scheme for large U, since the total angular momentum J of 
the ground-state multipole is purely composed of seven spins. 
In the limit of large A, all the coefficients converge to unity, 
since the ground state is essentially expressed by one elec- 
tron in the j=7/2 octet and the one-body multipole expression 
perfectly agrees with the Stevens operator in the one-electron 
case. It is interesting to mention that this fact is also related to 
Yb 3+ with one hole in the j=7/2 octet. 

It should be noted that around at X/U=0. 1^0.2, corre- 
sponding to the region for realistic values of actual materi- 
als, the coefficients for multipoles other than dipole rapidly 
increase, although the magnitude of the coefficient depends 
on individual values of U and A. In particular, it is impres- 
sive that the coefficient for quadrupole seems to be larger 
than we have naively expected, since the realistic situation 
with A/J7=0.1^0.2 has been believed to be well described by 
the LS coupling scheme. This result strongly suggests that 
quadrupole degree of freedom cannot be simply neglected in 



Gd compounds. The extent of the effect of quadrupole will be 
discussed in the next section. 

From a mathematical viewpoint, the relation eq. (13) just 
expresses the Wigner-Eckart theorem, but it contains physi- 
cally important message concerning the multipole. Since the 
J-multiplet of /"-electron system is given by the appropriate 
superposition of ri-body states, the Stevens operator is consid- 
ered to be given by the expression of multipole in the many- 
body form, in sharp contrast to the present definition of one- 
body form. However, for /"-electron systems, eq. (13) means 
that the expectation value of the multipole operator in the 
present definition is in proportion to the Stevens operator and 
the effect of interactions is included only in the proportional 
coefficient c. When we impose the orthonormal condition on 
the multipole operator, the difference in the normalization is 
absorbed in the proportional coefficient. Namely, the multi- 
pole in the present definition is equivalent to the Stevens op- 
erator, except for the proportional coefficient. Thus, it is pos- 
sible to use one of the definitions of multipole, one-body or 
many-body form, depending on the problem, but we point out 
some advantages to use the multipole in the one-body form. 

The Stevens operator is considered to be suitable for local- 
ized situation, which is well described by the LS coupling 
scheme. First we construct the Hund's-rule ground state char- 
acterized by spin moment S and angular momentum L for 
/"-electron systems. Then, we include the effect of spin-orbit 
coupling to form the multiplet characterized by total angular 
momentum J, which is given by J=\L—S\ and J=L+S for 
n<7 and n>7, respectively. In such a case, it seems to be nat- 
ural to define the multipole by the Stevens operator given by 
eq. (14). However, it is difficult to consider itinerant situation 
in the LS coupling scheme, since J is changed by charge fluc- 
tuations. In contrast, the present definition of multipole can be 
used both for itinerant and localized cases, since it is defined 
as the one-body density operator. This is one of advantages of 
the present definition of multipole. 

When we define multipole in the one-body density oper- 
ator, the effect of Coulomb interaction and spin-orbit cou- 
pling is automatically included in the proportional coefficient. 
Thus, in the present definition of multipoles, we can discuss 
the relative change of multipole moment in the wide range 
of Coulomb interaction and spin-orbit coupling from X/U=Q 
(LS coupling scheme) to X/U=oo (j-j coupling scheme). 
This is another advantage of the one-body definition of mul- 
tipole. It is true that from a symmetry viewpoint, we can con- 
tinue to use the Stevens operator even if we deviate from 
the LS coupling scheme, since the ground state multiplet is 
always characterized by the same total angular momentum 
when we change U and/or A. However, the coefficient c in 
eq. (13) is changed by the interactions, but this point is not 
included in the Stevens operator. This is easily understood 
when we recall the CEF potential for the case of n=2, in which 
the ground state multiplet is characterized by J=4. The CEF 
potentials are expressed by using the fourth- and sixth-order 
Stevens operators, but the sixth-order one vanishes in the j-j 
coupling regime, since the CEF term is originally given by the 
sum of one-electron electrostatic potential and the sixth-order 
terms do not appear for one electron in the j=5/2 states which 
are accommodated in the limit of large A. 

Here we should mention that there is a disadvantage con- 
cerning the limitation of the rank of multipole in the present 
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definition of the multipole in the one-body form. As men- 
tioned in the previous subsection, we can treat the multipoles 
up to rank 7 in the present multipole definition. In order to 
treat higher multipoles with fc>8, it is necessary to define the 
multipole in the many-body operator form, since the Stevens 
operator with the many-body form can treat all possible mul- 
tipoles up to rank 2 J in the multiplet characterized by J. Nev- 
ertheless we believe that the present multipole definition with 
one-body form keeps an advantage, since in actuality, it is 
very rare to treat higher-order multipoles with fc>8. 

3. Crystalline Electric Field States 

3.1 Effective Hamiltonian 

The local CEF state can be determined by the diagonaliza- 
tion of the local Hamiltonian eq. (1), but in order to under- 
stand how the CEF states are changed due to the competition 
between the Coulomb interaction and spin-orbit coupling, i.e., 
the LS and j-j coupling schemes, it is useful to derive the 
effective Hamiltonian H c g from eq. (1) with the use of the 
Stevens operator eq. (14). We express H c a as 

H cS (n, J) = B° 4 {n, J)(6° 4 + 50j) 



Bl(n,J){6l-210t) 
5 6 2 (n,J)(0 2 -6 6 6 ), 



(15) 



where Of, is the Stevens operator of which matrix elements 
for n and J are given by eq. (14) and B%(n, J) denotes the 
CEF parameter for the multiplet specified by n and J. 

The effects of U and A appear in the CEF parameters, 
B 4 (n, J), Bg(n, J), and Bq(ti, J), which are expressed by 

B° 4 (n,J) = k 4 (n,J)B° 4 , 

B°(n,J) = k 6 (n,J)Bl (16) 



Bi{n,J) = k 6 (n,J)B i 



6 ' 



where B 4 , Bg, and B\ for one / electron with £=3 are given 
by eq. (9) and k 4 (n, J) and fc 6 (n, J) are given by 



k 4 {n, J) = 4 n) /A, k 6 (n, J) = 7y*V 7/ , 



(17) 



respectively. Here pj and jj 1 ' are the so-called Stevens 
factors, which are coefficients appearing in the method of 
Stevens' operator equivalent. 29 Note that in general, pj and 

In) 

7} depend on the Coulomb interaction and spin-orbit cou- 
pling, since they are determined by the nature of the ground- 
state multiplet specified by J. 

For one / electron with 1=3, (3g and 7^ are given by 

2 4 



(») 



fa 



(18) 



11-45' '* 9 13 -33' 
respectively. For the ground state of n=l and J=7/2, in the 
limit of U=oo, i.e., in the LS coupling scheme, we obtain 
k 4 s (n, J)=kg S (n, J)=0, suggesting that the CEF potentials 
do not work for / states with J=S=1I2 in the LS cou- 
pling scheme. On the other hand, in the limit of A=oo, i.e., 
in the j-j coupling scheme, we obtain k? 4 ~ 3 (n, J)=3/7 and 
fcg ~ J (n, J)=l/7 from the Stevens factors 29 



P7/2 — 



15 • 77 



(7) 
T7/2 



13 ■ 33 • 63 



(19) 



Note that the absolute values of $?J 2 and Jy 2 in the j-j cou- 




Fig. 2. (Color online) Coefficients k^{n,J) and ka(n, J) for n=7 and 
J =7/2 obtained from the diagonalization of Hq + H so . 



pling scheme are equal to those for Yb 3+ ion with one hole in 
j=7/2 octet, although the signs are just inverted. 

For the intermediate coupling region in which both U and A 
are finite, we numerically evaluate kn(n, J) and kg(n, J) by 
deriving 7J c ff from the local Hamiltonian H\ oc . By using the 
ground state \n,J, M) of H so + Hq, we evaluate the matrix 
elements of the effective Hamiltonian as 



H eS {n, J; M, M') = (n, J, M\H CEF \n, J, M') 



(20) 



where we assume that \W\ is much smaller than both U and A. 
Since the matrix elements of the Stevens operator have been 
already listed for each value of J, we can numerically obtain 
Bl(n, J), Bg(n, J), and B%(n, J) from eq. (20) for any val- 
ues of U and A for a given value of n. 

In Fig. 2, we show k 4 (n, J) and ke(n,J) for n=7 and 
J =112 as functions of X/U. In order to draw this figure, 
we numerically evaluate the matrix element eq. (20) for 
|W|=10- 4 eV. Concerning U and A, for X/U < 1, we fix 
A as A=1.0 and increase U, while for X/U > 1, we fix U as 
[7=1.0 and increase A. As mentioned in Fig. 1, we have not 
found significant difference in the results, even when we indi- 
vidually change the values of U and A by keeping the ratio. 

As mentioned above, we obtain k 4 (n, J)=ke(n, J)=0 for 
small X/U, while we find k^n, J)=3fce(ri, J)=3/7 for large 
X/U. In the region of X/U=Q. 1^1.0, we observe rapid in- 
creases of both k 4 (n, J) and fce(ri, J) from zeros to the val- 
ues of the j-j coupling limit. In this sense, for actual values 
of A and U, we expect the appearance of the CEF effect even 
for the case of n=l . The quantitative argument of the CEF 
potentials will be discussed in the next subsection. 

Finally, we note the difference in the increasing behavior of 
k 4 (n, J) and k^n, J). Although k 4 (n, J) monotonically in- 
creases from zero to 3/7, kQ(n, J) does not change monoton- 
ically, but it forms a peak around at A/J7~1.0, In particular, 
k§(n, J) is slightly larger than 1/7, which is the value of the 
j-j coupling limit. In the region of X/U=0. 1^1.0, the effect 
of ke(n, J) is relatively larger than that of k^n, J). 

3.2 CEF parameters 

Before proceeding to the evaluation of quadrupole sus- 
ceptibility, it is necessary to determine the CEF parameters 
in -Hcef- For the purpose, as typical material, we consider 
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Fig. 3. CEF energy schemes for (a) n=2 and (b) n=7. As for the details of 
the parameters, see the main text. Here we note that we set W=— 0.69 meV 
and j/=1.08. 



filled skutterudite compounds. 4 As for the CEF parameters of 
PrOs4Sbi2, it has been confirmed that the ground state is 

-L(2) 

singlet state and the first CEF excited state is T 4 triplet with 
small excitation energy of about 10 K. 5-7 Note that since filled 
skutterudite compounds crystallize in the cubic structure with 
Th symmetry, two triplet states of and in Oh symmetry 
are mixed with each other, leading to two states. 10 

In order to reproduce the CEF energy scheme which has 
been experimentally determined, 5-7 we perform the optimiza- 
tion of W, x, and y in Hqef of eq. (1) for [7=1.0 eV. Note 
that A is set as the value which has been experimentally 
found. 32 For Pr 3+ ion, we use A=0.094 eV. After the opti- 
mization, we find 1^=— 0.69 meV, x=0.22, and y=1.08. We 
note that these CEF parameters determine the CEF poten- 
tials for one /-electron state of eq. (9), not for / 2 -states of 
J=4 of eqs. (15) and (16). In Fig. 3(a), we show the CEF 
energies vs. x for W=— 0.69 meV and y=1.08. The vertical 
broken line denotes the position of x=0.22, at which we find 
£( r 4 (2) ) - E(T+)=7.9K, E{r+ {1) ) - £(r+)=134.9K, and 
£(r+) -£(r+)=205.4K, where E(T^) denotes the eigenen- 
ergy of the CEF state characterized by the irreducible repre- 
sentation of l\y. These excitation energies reproduce well the 
experimental CEF energy scheme of PrOs4Sbi2. 5 ~ 7 

Here it is instructive to consider the CEF potentials in 
eq. (16) of the effective model for n=2 and J=4. After some 
algebraic calculations with the use of the above optimized val- 
ues of W, x, and y, we obtain i?2(n = 2, J = 4)=2.29 x 
10~ 2 K, B%{n = 2, J = 4)=1.33 x 10~ 3 K, and B%(n = 
2, J = 4)=1.38 x 10 _2 K. In the analysis of the experimen- 
tal results, 5 the CEF potentials for the states of n=2 and J=A 



have been found to be A 4 =2.37 x 10~ 2 K, ^ 6 =1.32 x 10~ 3 K, 
and Ag=1.08 x 10~ 2 K, where A 4 , A 6 , and A\ used in Ref. 5 
denote Bl{n = 2, J = 4), B%(n = 2, J = 4), and 
B\(n — 2, J = 4) in our notation, respectively. The CEF po- 
tentials for the states of n=2 and J=4 obtained by considering 
both finite spin-orbit coupling and Coulomb interactions are 
totally similar to those obtained in the LS coupling scheme 
in which the Coulomb interactions are assumed to be infi- 
nite, although there exist deviations between two groups of 
the CEF potentials due to the effect of finite Coulomb inter- 
actions. This fact indicates that the CEF states of n=2 are well 
described by the LS coupling scheme. 

Next we show the CEF energy scheme for Gd-based filled 
skutterudites by using the same CEF parameters. Since the 
CEF potentials work on / orbital, as easily understood from 
Hcef, there is no reason to change the CEF parameters when 
we replace the rare-earth ion with another one in the same 
environment due to the same ligand ions. Note that A is de- 
termined from the experimental value and we set A=0.18 eV 
for Gd ion. In Fig. 3(b), we show the CEF energies vs. T for 
n=7 with the use of the same parameters as those in Fig. 3(a), 
except for the spin-orbit coupling. Note that the magnitude of 
the CEF splitting is obviously smaller than that of Fig. 3(a), 
but we observe the CEF excitation in the order of a few 
Kelvin. In the present parameters, we find that at a;=0.22, the 
ground state is Tq 7 quartet and the excited states are two types 
of doublets. Note that two different doublets Tq and Tj 
in Oh symmetry are mixed with each other, leading to two 
doublet states. Note also that quartet in Oh symmetry is 
expressed by r^" 7 in Th symmetry. 

In the LS coupling scheme, there is no CEF splitting for the 
case of n=7, since the octet of J=7/2 is composed of 5=7/2 
and L=Q. However, as emphasized in the discussion of multi- 
pole operator, the actual situation is in the middle of the LS 
and j-j coupling schemes. Thus, in contrast to our knowledge 
on the basis of the LS coupling scheme, we observe the CEF 
splitting due to the effect of the j-j coupling component in 
the wavefunction. This point will be discussed again for the 
quadrupole susceptibility. 

4. Quadrupole Susceptibility 

Now we discuss the quadrupole susceptibility, which can 
be detected by the measurement of elastic constant C. In gen- 
eral, we consider the strain fields e u and e v which are coupled 



and Qvi=T- 2 j2_ y 2, respec 



with quadrupoles Q h 

tively. In the notation of the Stevens operator, these are ex- 
pressed by (2J| — J 2 — Jy)/2 and ^(Jj - J 2 )/2, respec- 
tively. 30 ' 31 

By following the standard formalism to calculate the elas- 
tic constant, we add the coupling term between strain and 
quadrupole to i?i oc , expressed as 



H = H 



loc 



QS 



H, 



QQ> 



(21) 



where Hq§ and Hqq denote quadrupole-strain interaction 
and inter-site quadrupole interaction terms, respectively. We 
express them as 



Hqs = -g^2{Qi 



(22) 
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and 



#QQ 



9 ^ ] {QuiQuj 
(id) 



f), 



(23) 



where 5 and 5', respectively, denote quadrupole-strain and 
inter-site quadrupole couplings and indicates the pair 

of nearest neighbor atomic sites. 

The elastic constant C is usually expressed with the use 
of quadrupole susceptibility \Q as C=—Ng 2 XQ, where N 
denotes the number of rare-earth ion in the unit volume. By 
considering only the strain field of e v , we estimate \Q m a 
simple mean-field approximation, in which we express -Hqq 
as 

H QQ, = -9' X! QviiQvj), (24) 
(id) 

where (• • ■ ) denotes the operation to take thermal average. 
Then, we obtain \Q m the mean-field approximation as 

X 



XQ = 



(25) 



1-9'X 

where \ lS the on-site quadrupole susceptibility determined 
from H\ oc , given by 

x- 7 E C V— V- IH[&-(&)]MI 2 - (26) 

Here we suppress the site dependence of Q v i, Z is the 
partition function given by Z=J2 n e~ En/T , E n denotes 
the eigenenergy of the rt-th eigenstate \n) of H\ oc , and 
(Q)=E n e- E " /T (n\Q\n)/Z. 

In Fig. 4(a), we depict — \q vs. T for n=2 correspond- 
ing to PrOs4Sbi2 for several values of g' with the use of 
[7=1.0 eV, A=0.094 eV (experimental value), and the CEF 
parameters determined in the previous subsection. For g'=0, 
we find a minimum in — xq at T«3K due to the suppres- 
sion of the Curie term. Then, at low enough temperatures, 
the quadrupole susceptibility is dominated by the Van Vleck 
term due to the off-diagonal term stemming from the quasi- 
degenerate Ti sing let and triplet, 33 - 34 which are charac- 
teristic to PrOs4Sbi2- When we increase the value of g 1 up to 
1 K as shown in Fig. 4(a), we do not find significant changes 
in the temperature dependence, although the absolute value is 
slightly suppressed. 

Here we remark that the characteristic behavior of ~xq 
for g'=0 has been obtained in the analysis of the experimen- 
tal results in the LS coupling scheme. 33 ' 34 The experimen- 
tal results on the elastic constant for PrOs4Sbi2 have been 
well explained by such calculation results. Here it is empha- 
sized that we correctly reproduce the previous results by using 
the definition of multipole in the one-electron density opera- 
tor form. Note, however, that the absolute value of — xq for 
g'=0 is different from that in the previous result, although the 
shape of the temperature dependence agrees quite well with 
the previous one. An important origin of such difference is 
due to the difference in the wavefunction. Another origin is 
due to the normalization of the multipole operator, which will 
be commented later. In any case, for the purpose to reproduce 
the experimental results with the calculation one, we adjust 
the value of g, which is just a fitting parameter, since it is 
not determined solely from the theoretical calculations in the 
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Fig. 4. (Color online) Quadrupole susceptibility vs. temperature T for (a) 
n=2 and (b) n=l. As for the details of the parameters, see the main text. Note 
here that the unit of g' is Kelvin. 



present framework. Thus, in this paper, we do not take care of 
the difference in the absolute values between the present xq 
and the previous results. However, it is meaningful to consider 
the relative difference between the results for n=2 and n=7 in 
the same calculation framework. 

In Fig. 4(b), we depict ~xq vs - T for n=l corresponding to 
Gd-based filled skutterudite compounds for several values of 
g' with the use of U=l.Q eV, A=0.18 eV (experimental value), 
and the same CEF parameters as those for PrOr4Sbi2. The 
results are similar to those for Ce compounds with Tg quar- 
tet ground state. This point is intuitively understood from the 
similarity between Ce 3+ ion with one electron in the j=5/2 
sextet and Gd 3+ ion with one electron in the j=7/2 octet in 
the limit of the j-j coupling scheme, when quartet ground 
state appears due to the CEF potentials. 

Note that as easily understood from the difference in the 
vertical scales between Figs. 4(a) and 4(b), the absolute value 
of the quadrupole susceptibility for n=l is larger than that 
for n=2, mainly due to the contribution of the Curie term. 
When we increase g 1 in the order of 1 K, the magnitude of 
the quadrupole susceptibility is totally suppressed, but it is 
still possible to observe the softening in the elastic constant, 
which is determined by the quadrupole susceptibility. For the 
case of PrOs4Sbi2, the fitting to the experimental results has 
been done for g'=0, 33 ' 34 but at the present stage, we cannot 
determine the magnitude of g 1 in order to discuss the absolute 
value of — xq of Gd-based filled skutterudites. An important 
message which we emphasize here is that it is possible to ob- 
serve the softening of elastic constant even in Gd-based filled 
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skutterudite compounds contrary to our knowledge on the ba- 
sis of the LS coupling scheme. 

5. Discussion and Summary 

We have evaluated the quadrupole susceptibility for realis- 
tic values of Coulomb interaction and spin-orbit coupling for 
Pr- and Gd-based filled skutterudite compounds. In particu- 
lar, we have pointed out a possibility that the softening of the 
elastic constant can be observed for Gd compounds. 

Concerning previous research on elastic properties of Gd 
compounds, we notice the measurement of elastic constant 
in rare-earth hexaborides. 35 No anomaly has been reported in 
GdB 6 due to the elastic constant measurement, except for the 
signal at a temperature of the antiferromagnetic phase transi- 
tion. In this sense, it may be difficult to observe actually the 
softening of elastic constant in Gd compounds with the cubic 
structure. However, we believe that it is worthwhile to per- 
form more precise ultrasonic experiments in Gd compounds, 
in particular, for Gd-based filled skutterudites. 

Here we also remark that Eu 2+ also contains seven / elec- 
trons. Thus, we expect the appearance of the softening of elas- 
tic constant in Eu compounds. 36 Note, however, that valence 
fluctuation between divalent and trivalent states is considered 
to be strong in europium. It may be also difficult to detect the 
softening of elastic constant in Eu compounds, but we hope 
that precise ultrasonic experiments in Eu compounds will be 
done in near future. 

In the evaluation of the quadrupole susceptibility, we have 
not considered the competition among different types of mul- 
tipoles. Since we have concentrated only on the measurement 
of the elastic constant, the quadrupole susceptibility has been 
evaluated. However, for Gd compounds, we have frequently 
observed an antiferromagnetic phase in the low-temperature 
region. When the Neel temperature is so high, it may be diffi- 
cult to detect the softening of the elastic constant in the para- 
magnetic phase. In such a case, it is interesting to consider the 
anisotropy of magnetization in the antiferromagnetic phase 
due to the effect of quadrupole degree of freedom. It is one 
of future problems to estimate the extent of the anisotropy of 
magnetization in Gd compounds. 

If we intend to show explicitly the occurrence of the anti- 
ferromagnetic phase, it is necessary to perform the optimiza- 
tion of multipole susceptibility. 23 ' 27 Then, we conclude that 
the multipole state with the largest eigenvalue of the suscep- 
tibility matrix is realized when we decrease the temperature. 
Such calculations have been actually performed for the seven- 
orbital Anderson model with the use of a numerical renormal- 
ization group technique. 21 For the case of Gd ion, it has been 
shown that the dipole state exhibits the maximum eigenvalue 
and the eigenstate of the second largest eigenvalue is charac- 
terized by quadrupole. As for the extension of the calculations 
to the periodic systems, 27 it is another future issue. 

In this paper, we have considered the coupling between 
quadrupole and the strain for the elastic constant. However, 
from the symmetry viewpoint, it is necessary to include also 
the effect of hexadecapoles. In fact, for the explanation of 
the temperature dependence of the elastic constant of PrMg 3 
with T 3 doublet ground state, a crucial role of hexadecapole 
has been emphasized. 37 Also in Gd-based compounds, such 
higher-order multipoles may play some roles for the explana- 
tion of the behavior of elastic constant. If we will successfully 



observe the softening of the elastic constant in Gd-based filled 
skutterudites, it may be necessary to include the effect of hex- 
adecapole to explain precisely the temperature dependence of 
the elastic constant. It is also another future problem. 

Finally, we provide a comment on the normalization of the 
multipole operator. In the present paper, since we have con- 
sidered only the quadrupole susceptibility, it is not necessary 
to pay attention to the normalization of the multipole opera- 
tor. However, if we consider the diagonalization of the multi- 
pole susceptibility matrix, it is inevitable to define the normal- 
ization. A natural way is to redefine each multipole operator 

~ (k) ~ (k') i 

so as to satisfy the orthonormal condition of Tr{T 7 'T\ ' } 
=<5fcfc"577'- 38 Note again that we have not considered such re- 
definition in the present paper. 

In summary, we have found that the j-j coupling compo- 
nent significantly appears in the ground-state wavefunction 
for n=l with J=7/2 even for realistic values of Coulomb in- 
teraction and spin-orbit coupling. With the use of the CEF pa- 
rameters determined from PrOs4Sbi2, we have evaluated the 
quadrupole susceptibility for Gd-based filled skutterudites. 
We have emphasized that the softening of the elastic constant 
can be detected in Gd compounds. We expect that ultrasonic 
experiments will be performed in Gd-based filled skutterudite 
compounds in a future. 
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